Cancer Regression Analysis

Data Loading and Initial Exploration

library(tidyr)
library(lmtest)
library(MASS)
library(car)
library(dplyr)
library(ggplot2)
library(glmnet)
library(gridExtra)

cancer <- read.csv("/Users/kylie/Documents/Math372/cancer_reg.csv")
death.rate <- cancer$target_deathrate

hist(death.rate, main = "Death Rate Distribution", 
     xlab = "Death Rate", col = "lightblue", border = "white")

Distribution of Cancer Death Rates

Data Preprocessing

We removed redundant and highly correlated features, then separated geographic information into county and state variables.

# Remove unwanted features
features_to_remove <- c("pctsomecol18_24", "pctprivatecoveragealone",
                        "pctpubliccoveragealone", "pctemployed16_over", "binnedinc",
                        "avganncount", "avgdeathsperyear", "popest2015")
cancer_clean <- cancer[, !(names(cancer) %in% features_to_remove)]

# Separate geography
cancer_clean <- cancer_clean %>%
  separate(geography, into = c("county", "state"), sep = ", ", remove = FALSE)

features_to_remove <- c("geography", "county")

# Convert state to factor
cancer_clean$state <- as.factor(cancer_clean$state)
cancer$state <- cancer_clean$state
cancer_clean <- cancer_clean[, !(names(cancer_clean) %in% features_to_remove)]

num_states <- length(levels(cancer_clean$state))

The dataset contains 51 states and 3047 counties.

Variable Selection with LASSO

cancer.rate <- c("incidencerate")
cancer_clean_death_pred <- cancer_clean[, !(names(cancer_clean) %in% cancer.rate)]
death_pred_matrix <- model.matrix(target_deathrate ~ ., data = cancer_clean_death_pred)

set.seed(1)
grid.lambda <- 10^seq(10, -2, length = 100)
x <- death_pred_matrix
y <- cancer_clean$target_deathrate
train <- sample(1:nrow(x), round(0.9 * nrow(x)))
test <- (-train)

x.train <- x[train, ]
x.test  <- x[test, ]
y.train <- y[train]
y.test  <- y[test]

cv.lasso <- cv.glmnet(x.train, y.train, alpha = 1)
best.lambda.lasso <- cv.lasso$lambda.min
lasso.pred <- predict(cv.lasso, s = best.lambda.lasso, newx = x.test)
mspe.lasso <- mean((lasso.pred - y.test)^2)

plot(cv.lasso)

LASSO Cross-Validation Results

The LASSO model achieved an MSPE of 511.9.

# Get LASSO coefficients at best lambda
lasso_coefs <- coef(cv.lasso, s = best.lambda.lasso)
selected_var_names <- rownames(lasso_coefs)[which(lasso_coefs != 0)]
selected_var_names <- selected_var_names[selected_var_names != "(Intercept)"]

# Create data frame with selected variables
x_no_intercept <- x[, -1]
available_vars <- colnames(x_no_intercept)
valid_selected_vars <- selected_var_names[selected_var_names %in% available_vars]

selected_data <- data.frame(
  target_deathrate = y,
  x_no_intercept[, valid_selected_vars]
)

# Fit OLS on selected variables
model <- lm(target_deathrate ~ ., data = selected_data)

Model Diagnostics

par(mfrow = c(2, 2))
plot(model)

Initial Model Diagnostics

par(mfrow = c(1, 1))

residuals <- model$residuals
sd_residuals <- sd(model$residuals)
mean_residuals <- mean(model$residuals)
stand_residuals <- (residuals - mean_residuals) / sd_residuals
shapiro_result <- shapiro.test(model$residuals)
bp_results <- bptest(model)

The Shapiro-Wilk normality test rejected the null hypothesis therefore our normality assumptions have not been met. Similarly, Breusch-Pagan constant variance test also rejected the null hypothesis therefore our variance was found to be heteroscedastic. We found heavy tails on the QQ plot therefore decided to subset our data into three subsets (left, center, right).

Data Subsetting and Transformation

model_center_subset <- lm(target_deathrate ~ ., 
                          data = selected_data[abs(stand_residuals) < 2, ])

Our center subset was still non normal and had a non constant variance so we decided to apply a Box-Cox transformation.

Box-Cox Transformation

bc_result <- boxcox(model_center_subset, plotit = TRUE)

Box-Cox Transformation Analysis

optimal_lambda <- bc_result$x[which.max(bc_result$y)]
max_loglik <- max(bc_result$y)
ci_threshold <- max_loglik - qchisq(0.95, 1)/2
ci_indices <- which(bc_result$y >= ci_threshold)
lambda_ci <- range(bc_result$x[ci_indices])

Our optimal lambda is 0.707 with a 95% CI for lambda: [0.545, 0.869]. Therefore, we will apply a square root transformation. We will also apply square root transformations to our left and right subset, but we will focus our analysis on our center subset because it has the bulk of our data, >2800 observations.

Transformed Models

# MAIN MODEL --> where most of our data is >2800 data points focus analysis on this model
model_center_subset_sqrt <- lm(sqrt(target_deathrate) ~ ., 
                               data = selected_data[abs(stand_residuals) < 2, ])

# Left subset
model_left_subset <- lm(sqrt(target_deathrate) ~ ., 
                        data = selected_data[stand_residuals < -2, ])

# Right subset
model_right_subset <- lm(sqrt(target_deathrate) ~ ., 
                         data = selected_data[stand_residuals > 2, ])

par(mfrow = c(2, 2))
plot(model_center_subset_sqrt)

Transformed Center Model Diagnostics

par(mfrow = c(1, 1))

Research Question 1: Insurance Coverage Effects

Does the Percent of Private/Public Healthcare Coverage effect cancer related deaths in those counties?

main_coef <- coef(model_center_subset_sqrt)[c("pctempprivcoverage", "pctpubliccoverage")]
left_coef <- coef(model_left_subset)[c("pctempprivcoverage", "pctpubliccoverage")]
right_coef <- coef(model_right_subset)[c("pctempprivcoverage", "pctpubliccoverage")]

knitr::kable(
  data.frame(
    Subset = c("Center", "Left", "Right"),
    Private = c(main_coef[1], left_coef[1], right_coef[1]),
    Public = c(main_coef[2], left_coef[2], right_coef[2])
  ),
  caption = "Insurance Coverage Coefficients by Subset",
  digits = 5
)

Insurance Coverage Coefficients by Subset

Subset	Private	Public
Center	0.01467	0.03563
Left	0.02073	0.06402
Right	0.01842	0.03829

In our center subset, the coefficient estimate for percent private employer-based coverage is 0.01467, with a p-value of 0.0356, indicating that this predictor is statistically significant. This means that for each 1% increase in private employer-based insurance coverage, the square root of the cancer death rate increases by 0.0147 units, holding all other variables constant.

The coefficient estimate for percent public insurance coverage is 0.03563, with a p-value < 2e-16, making it highly significant. This positive association indicates that counties with higher public insurance enrollment tend to have higher cancer death rates. Because public coverage (Medicaid, Medicare) is most common in low-income, older, and medically vulnerable populations, this likely reflects underlying socioeconomic and health-status differences rather than any effect of insurance type itself. For each 1% increase in public insurance coverage, the square root of the cancer death rate increases by 0.0356 units, holding all else constant.

Across all three subsets (center, left, and right), public coverage consistently has a larger coefficient than private coverage, implying that counties with higher proportions of public insurance tend to have higher observed cancer mortality. These coefficients likely reflect population-level risk factors—such as age distribution, baseline health, and socioeconomic disadvantage—rather than the quality of the insurance plans themselves.

ggplot(cancer_clean, aes(x = pctempprivcoverage, y = target_deathrate)) +
  geom_point(aes(color = pctpubliccoverage), alpha = 0.7, size = 2) +
  scale_color_gradient(low = "pink", high = "lightblue") +
  geom_smooth(method = "lm", color = "darkgreen", se = TRUE) +
  labs(
    title = "Private and Public Coverage vs Cancer Death Rate",
    x = "Private Insurance Coverage (%)",
    y = "Cancer Death Rate (per 100,000)",
    color = "Public Coverage (%)"
  ) +
  theme_minimal(base_size = 14)

Private and Public Coverage vs Cancer Death Rate

This graphic compares private and public insurance coverage vs cancer death rate. The x-axis shows the percentage of people in a county with private insurance, and the y-axis shows the number of cancer deaths per 100,000 in a county. Each point represents a county, and the color of each point represents whether the county has high (light blue) or low (pink) public coverage. Using this, we are able to see how counties are clustered based upon the combination of public and private coverage. The green line shows the linear regression line of private coverage vs death rate. This shows a downward slope, meaning higher private coverage is associated with lower death rate. Notice the counties with low public coverage have lower death rates. Blue points, counties with high public coverage, cluster higher on the y-axis, confirming higher death rates. This signals that it is likely that public coverage is more common in high-risk and/or underserved populations.

cancer_clean <- cancer_clean %>%
  mutate(
    private_group = cut(
      pctprivatecoverage,
      breaks = quantile(pctprivatecoverage, probs = c(0, .33, .66, 1), na.rm = TRUE),
      include.lowest = TRUE,
      labels = c("Low", "Medium", "High")
    ),
    public_group = cut(
      pctpubliccoverage,
      breaks = quantile(pctpubliccoverage, probs = c(0, .33, .66, 1), na.rm = TRUE),
      include.lowest = TRUE,
      labels = c("Low", "Medium", "High")
    )
  )

p1 <- ggplot(cancer_clean, aes(x = private_group, y = target_deathrate)) +
  geom_boxplot(fill = "lightpink") +
  labs(
    title = "Cancer Death by Priv Insurance",
    x = "Private Coverage Group",
    y = "Cancer Death Rate (per 100,000)"
  ) +
  theme_minimal(base_size = 14)

p2 <- ggplot(cancer_clean, aes(x = public_group, y = target_deathrate)) +
  geom_boxplot(fill = "lightblue") +
  labs(
    title = "Cancer Death by Pub Insurance",
    x = "Public Coverage Group",
    y = "Cancer Death Rate (per 100,000)"
  ) +
  theme_minimal(base_size = 14)

grid.arrange(p1, p2, ncol = 2)

Cancer Death Rate by Insurance Coverage Level

Here we can see the real trends in the data. Beginning with cancer death rate with private insurance coverage. As private coverage increases from low to medium to high, the median cancer death rate appears to decrease. The interquartile range is somewhat consistent across all the groups, but the low private coverage group has counties with high death rates. Counties with higher private insurance coverage tend to have lower cancer death rates, suggesting that private coverage can lead to better access to preventative care, early detection, and treatment. However, notice there are still outliers, indicating there may be other factors involved such as income, urbanization, and healthcare infrastructure. As public coverage increases from low to medium to high, the median cancer death rate increases. The interquartile range is wider for the public insurance groups, with several counties having much higher death rates. We can infer that high public coverage is associated with higher cancer death rates; however, this does not mean public coverage causes death, rather this coverage is in high-risk counties. Overall, we can see that these are inversely related in terms of cancer death rate. We can conclude that insurance coverage type is correlated with underlying county characteristics that affect mortality, rather than directly influencing cancer death rates.

Research Question 2: Education Effects

Does the Education background of County citizens effect cancer related deaths in those Counties?

Looking at the coefficients, there appears to be an association with highest education level earned.

vars <- cancer_clean[, c("target_deathrate", "pcths25_over", "pctbachdeg25_over", 
                         "pcths18_24", "pctbachdeg18_24", "pctnohs18_24")]
colnames(vars) <- c("Death Rate", "HS (25+)", "Bach (25+)", 
                    "HS (18–24)", "Bach (18–24)", "No HS (18–24)")
pairs(vars, pch = 1, main = "Scatterplot Matrix of Death Rate vs Education Levels by Age")

Scatterplot Matrix of Death Rate vs Education Levels by Age

p1 <- ggplot(cancer_clean, aes(pctbachdeg25_over, target_deathrate)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", color = "darkgreen") +
  labs(title = "Population 25+ with Bachelor's Degree vs Death Rate",
       x = "Percent of Population", y = "Death Rate") +
  theme_minimal()

p2 <- ggplot(cancer_clean, aes(pctbachdeg18_24, target_deathrate)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", color = "darkgreen") +
  labs(title = "Percent of Population 18-25 with Bachelor's Degree vs Death Rate",
       x = "Percent of Population", y = "Death Rate") +
  theme_minimal()

p3 <- ggplot(cancer_clean, aes(pcths18_24, target_deathrate)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", color = "darkgreen") +
  labs(title = "Percent of Population 18-24 with Highschool Diploma vs Death Rate",
       x = "Percent of Population", y = "Death Rate") +
  theme_minimal()

p4 <- ggplot(cancer_clean, aes(pctnohs18_24, target_deathrate)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", color = "darkgreen") +
  labs(title = "Percent of Population 18-24 with No Highschool Diploma vs Death Rate",
       x = "Percent of Population", y = "Death Rate") +
  theme_minimal()

p5 <- ggplot(cancer_clean, aes(pcths25_over, target_deathrate)) +
  geom_point(alpha = 0.4) +
  geom_smooth(method = "lm", color = "darkgreen") +
  labs(title = "Percent of Population 25+ with Highschool Diploma vs Death Rate",
       x = "Percent of Population", y = "Death Rate") +
  theme_minimal()

grid.arrange(p1, p2, p3, p4, p5, ncol = 2)

Education Attainment vs Cancer Death Rate

As the percentage of the population 25 years or older with the highest level of education obtained is Bachelor’s degree increases within a county, the mean per capita cancer mortality rate decreases. As the percentage of the population 25 years or older with the highest level of education achieved is a High school degree increases within a county, the mean per capita cancer mortality rate increases. This could be due to factors such as access to healthcare coverage or accessibility to information regarding healthier lifestyles and cancer detection measures.

main_ed <- coef(model_center_subset_sqrt)[c("pctnohs18_24", "pctbachdeg18_24", 
                                             "pcths18_24", "pcths25_over", "pctbachdeg25_over")]
left_ed <- coef(model_left_subset)[c("pctnohs18_24", "pctbachdeg18_24", 
                                     "pcths18_24", "pcths25_over", "pctbachdeg25_over")]
right_ed <- coef(model_right_subset)[c("pctnohs18_24", "pctbachdeg18_24", 
                                       "pcths18_24", "pcths25_over", "pctbachdeg25_over")]

knitr::kable(
  data.frame(
    Variable = c("No HS (18-24)", "Bach (18-24)", "HS (18-24)", "HS (25+)", "Bach (25+)"),
    Center = main_ed,
    Left = left_ed,
    Right = right_ed
  ),
  caption = "Education Coefficients by Subset",
  digits = 5,
  row.names = FALSE
)

Education Coefficients by Subset

Variable	Center	Left	Right
No HS (18-24)	-0.00143	-0.02678	-0.01374
Bach (18-24)	-0.01020	-0.04733	-0.00711
HS (18-24)	0.00612	0.00545	0.01451
HS (25+)	0.00964	0.01465	0.00799
Bach (25+)	-0.04821	-0.01421	0.01060

Comparing the coefficients directly, we check to see if they logically make sense to interpret. It makes sense that as the percent of the population with only a high school degree increases, death rate is expected to increase as well. Similarly, while the percent with bachelor’s degree increases, death rate is expected to decrease. Less intuitively, as the percent with no high school diploma 18-24 increases, death rate is expected to decrease, but this most likely is due to this group being younger (risk for cancer increases as you age), and the fact that this feature’s data is so close together. Interestingly, for ages 25 and over, the coefficients are what we would expect with high school only correlated with increasing death rate and a bachelor’s degree correlated with decreasing death rate, but for our right subsection of the data, as the percent of those with a bachelor’s degree 25 and over increases, death rate increases.

Research Question 3: Geographic Effects

How does geographical location (state) affect cancer related death rates in counties?

# Define regions
RegionFW <- c("California", "Nevada", "Oregon", "Washington", "Alaska", "Hawaii")
RegionRM <- c("Idaho", "Utah", "Colorado", "Wyoming", "Montana")
RegionSW <- c("Arizona", "New Mexico", "Texas", "Oklahoma")
RegionPl <- c("North Dakota", "South Dakota", "Nebraska", "Kansas", "Minnesota", "Iowa", "Missouri")
RegionGL <- c("Wisconsin", "Illinois", "Indiana", "Michigan", "Ohio")
RegionSE <- c("Louisiana", "Arkansas", "Mississippi", "Alabama", "Georgia",
              "Florida", "Tennessee", "South Carolina", "North Carolina",
              "Kentucky", "West Virginia", "Virginia")
RegionME <- c("Maryland", "District of Columbia", "Delaware", "Pennsylvania",
              "New Jersey", "New York")
RegionNE <- c("Connecticut", "Rhode Island", "Massachusetts", "New Hampshire", "Vermont", "Maine")

# Make region column
cancer_clean$states <- as.character(cancer$state)
cancer_clean$region <- NA

# Get region for each row
for (i in 1:nrow(cancer_clean)) {
  if (cancer_clean$states[i] %in% RegionNE) {
    cancer_clean$region[i] <- "New England"
  } else if (cancer_clean$states[i] %in% RegionME) {
    cancer_clean$region[i] <- "MidEast"
  } else if (cancer_clean$states[i] %in% RegionSE) {
    cancer_clean$region[i] <- "South East"
  } else if (cancer_clean$states[i] %in% RegionGL) {
    cancer_clean$region[i] <- "Great Lakes"
  } else if (cancer_clean$states[i] %in% RegionPl) {
    cancer_clean$region[i] <- "Plains"
  } else if (cancer_clean$states[i] %in% RegionSW) {
    cancer_clean$region[i] <- "South West"
  } else if (cancer_clean$states[i] %in% RegionRM) {
    cancer_clean$region[i] <- "Rocky Mountains"
  } else if (cancer_clean$states[i] %in% RegionFW) {
    cancer_clean$region[i] <- "Far West"
  }
}

# Check for any unassigned regions
unassigned <- cancer_clean %>% filter(is.na(region)) %>% distinct(states)
if (nrow(unassigned) > 0) {
  cat("Warning: The following states were not assigned to a region:\n")
  print(unassigned$states)
}

# Remove rows with missing regions (if any)
cancer_clean <- cancer_clean %>% filter(!is.na(region))
cancer_clean$region <- as.factor(cancer_clean$region)

ggplot(cancer_clean, aes(x = region, y = target_deathrate, fill = region)) +
  geom_boxplot() +
  labs(title = "Cancer Death by Region",
       x = NULL,
       y = "Cancer Death Rate (per 100,000)") +
  theme_minimal(base_size = 12) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1),
        legend.position = "none")

Cancer Death Rate by Region

These box plots seem to support our hypothesis that there is a difference between regions in terms of cancer deaths. The Rocky Mountains, Far West and New England have noticeably lower median death rates and relatively smaller interquartile ranges, which implies there isn’t much disparity between counties in their states. The South East and parts of the South West have some of the highest median death rates, with wider IQRs and several very high outliers, this could suggest that there are more disparities in death rates between counties within these states.

ggplot(cancer_clean, aes(x = state, y = target_deathrate, fill = region)) +
  geom_boxplot() +
  labs(title = "Cancer Death by State",
       x = NULL,
       y = "Cancer Death Rate (per 100,000)") +
  facet_grid(cols = vars(region), scales = "free_x") +
  theme_minimal(base_size = 10) +
  theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

Cancer Death Rate by State within Region

library(usmap)

# Get population back as a variable so we can aggregate the death rate correctly
cancer_clean$popest2015 <- cancer$popest2015

state_death <- cancer_clean %>%
  group_by(state) %>%
  summarise(
    # Keep region for later visual
    region = first(region),
    # Sum all county populations
    total_population = sum(popest2015),
    # Calculate the death rate for the entire state
    weighted_deathrate = sum(target_deathrate * popest2015) / sum(popest2015)
  )

# Repeat for region
region_death <- cancer_clean %>%
  group_by(region) %>%
  summarise(
    total_population = sum(popest2015),
    weighted_deathrate = sum(target_deathrate * popest2015) / sum(popest2015)
  )

The map reinforces these ideas, states in the South East and Midwest (more specifically parts of Appalachia and the deep South) seem to have much higher cancer mortality. We also see a kind of split in the South West between quite low mortality in states like New Mexico and Arizona, alongside one of the highest mortality rates for the whole country, Oklahoma. Like the boxplots, the Rocky Mountains, far West and part of New England have consistently lower mortality rates.

# Fix for the legend.justification.inside error
library(usmap)

# Create the base map
p <- plot_usmap(data = state_death, values = "weighted_deathrate") +
  scale_fill_gradientn(
    name = "Death Rate\n(per 100,000)",
    colors = topo.colors(10)
  )

# Add theme without problematic elements
p <- p + theme_void() +
  theme(
    legend.position = "right",
    legend.title = element_text(size = 10),
    legend.text = element_text(size = 8),
    plot.margin = margin(10, 10, 10, 10)
  )

# Print the plot
p

Cancer Death Rate by State (US Map)

These plots demonstrate there is substantial geographic variation in cancer mortality. These differences could point to underlying patterns in socioeconomic status as well as behavioral risk factors like smoking and environmental conditions. Further areas of exploration could be narrowing in on states with high variability by county or high overall cancer mortality and completing further analysis at the county level.

Regression: Because our model used dummy variables, the coefficient of each state represents the expected difference in cancer mortality rate when going from the reference variable, Alabama to that state, holding all other variables constant. Many of the state variables are significant to the model, showing that geography is a useful predictor in cancer death rates. The states with the greatest negative coefficients include Hawaii (-1.527), Utah(-1.169), Arizona (-1.354), Colorado(-0.996) and New Mexico (-0.944), it is notable that after the top 5, the next most negative coefficient is considerably less negative with California having a coefficient of -0.628. These states have considerably lower cancer mortality rates compared to Alabama. The states with the most significant and greatest positive coefficients include: Kentucky (0.660), Oklahoma (0.528), Tennessee (0.369), Arkansas (0.324) and Louisiana (0.309). These states have higher than expected cancer rates, even when compared to Alabama, which in our exploratory analysis seemed to have a relatively high death rate. While these coefficients could point to environmental conditions having an effect on cancer death rates, they also may point to larger issues concerning poverty alongside healthcare outreach and access rather than being a direct influence of the physical location.